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Abstract 

We compute the bulk viscosity of a gas of pions at temperatures below the QCD crossover tem- 
perature, for the physical value of m^, to lowest order in chiral perturbation theory. Bulk viscosity 
is controlled by number-changing processes which become exponentially slow at low tempera- 
tures when the pions become exponentially dilute, leading to an exponentially large bulk viscosity 
C ~ {F§/ml) exp(2m^/r), where Fq ~ 93MeV is the pion decay constant. 
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I. INTRODUCTION 



One of the most prominent discoveries of the heavy ion program at RHIC has been the 
success of hydrodynamics [1] with a zero [2] or very small [3] viscosity. Though the exact 
value of the viscosity cannot yet be extracted due to uncertainties in the initial state and 
other effects, it is a robust result that the viscosity near the QCD crossover temperature is 
small, 77/s < 0.5 [3]. On the other hand, perturbative calculations show that the viscosity 
to entropy ratio rj/s at high temperatures T ^ 1 GeV, where perturbation theory should 
work, is significantly higher [4]. Both theoretical [5] and data-driven [6] analyses of the pion 
gas indicate that rj/s also rises at low temperatures, suggesting that the relative viscosity 
bottoms out near the crossover [7], similar to the behavior in conventional fluids [8]. 

The bulk viscosity is also expected to be important in the hydrodynamics of heavy ion 
collisions [9]. Bulk viscosity vanishes for a conformal system, a good approximation to 
QCD at high temperatures; therefore the bulk viscosity to entropy ratio C/-^ is small at 
high temperatures [10]. Near the crossover temperature QCD is very far from conformal, 
as indicated by the peak in (e — 3P)/T^ [11, 12], and it is expected that (/s may display 
a peak at this scale [13, 14]. At lower temperatures QCD is well described by a pion gas. 
Existing studies of pion gases indicate that the bulk viscosity falls away at low temperatures 
[5, 6]. This suggests that the ratio (/s shows the opposite behavior of r^/s, peaking near the 
transition and falling ofi^ to either side [15]. 

However, previous analyses of the bulk viscosity of a pion gas have been very incom- 
plete. In particular, neither standard reference [5, 6] considers number changing processes. 
But such processes are essential to the relaxation of particle number to equilibrium and 
frequently control the bulk viscosity, as emphasized by Jeon [16]. Therefore we believe that 
what the calculations in the literature we computing was not really the bulk viscosity of a 
pion gas, but the constant for a relaxation process which treated kinetic but not chemical 
equilibration. To make a fair comparison with the calculations of C/s at higher temperatures 
one should compute the true bulk viscosity of a pion gas at low temperatures. When the 
bulk viscosity calculated in this way becomes large, it indicates that the pion gas will lose 
chemical equilibrium, a physically interesting property. 

In this paper we will provide a calculation of the bulk viscosity of a pion gas, including the 
relaxation via number changing reactions to chemical equilibrium. We will work to lowest 
nontrivial order in chiral perturbation theory, the effective theory of low energy pions. That 
is, we will make an expansion to lowest order in m7r/47rFo and T/AttFq, treating T/ttIt^ as a 
free parameter of order 1. (Here Fq is the pion decay constant, rriTr is the pion mass treating 
the ttq and tt^ as degenerate, and T is the temperature as usual.) Our treatment is therefore 
only valid at temperature scales low enough that there are almost no resonances (such as p 
mesons) and few kaons relative to pions; we will not try to extrapolate close to the crossover 
temperature. 

In the next section of the paper we will review the physics of bulk viscosity in a gas 
of relativistic, massive, weakly coupled bosons, emphasizing the role played by number 
changing processes. We show that the bulk viscosity is controlled by rrij^/T and by the 
rate of number changing processes. In Section 111 we present the calculation of the number 
changing rate within chiral perturbation theory. Our numerical results and conclusions 
are presented in Section IV, but can be summarized here. We flnd that, as temperature 
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falls, number changing processes become less efficient and the bulk viscosity actually grows, 
scaling as C/s ~ FqT'^^ for T ~ and C/-^ ~ F^T^^^'^niTr^^^'^ cxp{3mT^/T) for ^ T. 
Therefore the behavior of bulk viscosity is not the opposite of the behavior of shear viscosity, 
and in particular both the bulk viscosity to entropy ratio and the bulk viscosity itself diverge 
exponentially in the low temperature limit. 



II. KINETIC DESCRIPTION OF BULK VISCOSITY 



By definition, bulk viscosity C is a reduction of the pressure in an expanding system, and 
increase in pressure in a contracting system, proportional to the rate of volume change,^ 

P = Pe,-CV-t; = Pe,-C^- (2.1) 

This arises because the volume change induces a departure from equilibrium, which in turn 
modifies the pressure. To see how this occurs for a pion gas, we need to describe the 
system in terms of a calculable approximation scheme. Since physical QCD is near the 
chiral limit, the pion is a pseudo-Goldstone boson of the (spontaneously but also explicitly 
broken) chiral symmetry, and pions are therefore weakly coupled at low momenta and well 
described by chiral perturbation theory (see for instance [17, 18]). Weak coupling means 
that thermal pions will have well defined quasiparticles which will be well described by 
Boltzmann equations. Defining the species sum and integration 

the pressure is related to the occupancy of species a at momentum p, fa{p), via 

\ [ 2//a(p) . (2.3) 

Jap 

faip) in turn evolves according to the Boltzmann equation 

2E^^I^ + 2p . = -C[f] = -Ce,a«tic[/] - Cinel[/] , (2.4) 

with C[f] the collision operator, which we discuss in more detail below. 

The lefthand side of the Boltzmann equation drives the system from equilibrium. Since 
the bulk viscosity involves one spacetime gradient, we can find it by expanding the Boltz- 
mann equation to first order in gradients; since the lefthand side is explicitly first order in 
gradients, we may substitute fa{p,t) with its equilibrium form 



/o = exp 



Ep-v-p 

T 



1 ■ (2.5) 



^ When we write noncovariantly we implicitly work in the instantaneous local rest frame. We use boldface 
p, V for vectors and normal letters p, v for their magnitudes; P is always the pressure. 
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We take the energy to be = •y/p^~+mf, meaning that we will neglect interaction self- 
energy corrections in comparison to the explicit pion mass. Clearly this treatment does 
not allow lis to consider QCD in the strict chiral symmetry limit, where interaction effects 
are the only thing which lead to modified dispersion. It would be interesting to return to 
this case in the future, but we expect it to be rather subtle; for instance, the lowest order 
interaction effect actually does not change the dispersion relation [19, 20], and the next order 
only shifts the speed of propagation away from the speed of light [19, 21], which we believe 
also does not lead to nonvanishing bulk viscosity. Therefore interaction effects appear to 
arise at a rather high order in (T/47rFo)^. Therefore interaction effects can be neglected for 
T < m^, which is what we are considering. In this case, explicitly evaluating the lefthand 
side of the Boltzmann equation yields 

^Ep^^ + 2p ■ - 2/o(l+/o) (^^^ + -Jr^^VJj . (2.6) 

We are interested in the case diVj = V ■ v. The temperature changes with time because 
expansion causes cooling; at first order in gradients the time dependence of the temperature 
has its usual equilibrium relation to V • v, dT/dt — —ciT'V • v with = dP/de the squared 
speed of sound [10]. Therefore the lefthand side of the Boltzmann equation is 

2/o(l+/o) ^ V-v. (2.7) 

This "source" for departure from equilibrium has no net energy content. To see this, note 
that 

j %h{p) , 6 = / 2Elfoip) , (2.8) 
dP _dP/dT _ /„p¥#/o(l+/o) 
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de de/dT /^^2£;2^/o(l+/o) 



(2.9) 



and therefore 



/ Efo{l+fo)% = cl [ i?/o(l+/o)^, (2.10) 

J ap '^-^ J ap ^ 

which shows that there is no energy content for the "source" for departure from equilibrium. 
That this occurs is just a check that we have correctly identified the time dependence of the 
temperature. However, the "source" does carry a net particle number, namely 



dn 
'dt 



»2 _ Sc^E^ 



LHS 'Jap 



^ /o(l+/o)2 " 3^:^" V • ^ . (2.11) 



This means that expansion leaves excess pions, relative to the equilibrium number at the 
given energy density. The relaxation of this excess particle number controls equilibration 
and bulk viscosity. 

Next we turn to the collision term. At lowest (fourth) order in , the collision term 
contains only elastic tttt -H- tttt scattering. Such terms drive fa{p) towards its equilibrium 
form except that they cannot change total particle number. That is, there is no solution to 
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the linearized Boltzmann equation with Eq. (2.7) on the Icfthand side and only tttt f-)- tttt 
collision processes on the righthand side, since the lefthand side includes a change to the 
net particle number while the righthand side cannot change particle number. Therefore a 
calculation involving only number-conserving processes is incomplete and inconsistent, as 
emphasized by Jeon [16] in the context of scalar theory.^ Therefore we must include 
as well the lowest order number-changing process. Since QCD is parity symmetric but the 
pion is a parity-odd scalar, all interaction terms are even in the pion field and the lowest 
order kinematically allowed number- changing process is tttt -H- TTTTTTTr. 

At this point there is a simplification. As in the case of scalar theory [16] (but unlike 
the case of weakly coupled QCD [10]), number-changing processes are much less efficient 
than number-conserving processes in a pion gas. Number conserving processes drive the 
nonequilibrium distribution function f{p) = fo + 5f towards an almost-equilibrium form, 
but with a chemical potential for particle number, 

/».(exp ^-;;-;- -l)"'. (2.12) 

Here 5T is determined by the condition that the energy content of is the same as the 
energy content of /q. But number conserving processes cannot lead to the relaxation of /i 
towards zero, because the elastic collision term vanishes if /(p) = fn{p)- 

^' Jbp',ckdk' 

X - [/ ^ (1+/)] ) (2.13) 

as the gain term oc f{p) and the loss term oc (l-|-/(p)) cancel. Therefore f{p) will equal 
Pl^s ^ small correction. The value of /i will dominate the pressure shift. 
We cannot make the substitution f{p) = f^{p) in the elastic part of the collision operator. 
But if we consider the integral j^^ of Eq. (2.4), then the integral over Cdastic exactly vanishes, 
independent of the form of f{p). We can approximate f{p) = fnip) in the smaller inelastic 
part of C, yielding 

V V f fo{l+fo / ~ = / (-Cinel[/1) = -a„el ■ (2-14) 

J ap '^-^ J ap 

There are two contributions to this collision term. One contribution arises when p = pi is 
one of the four pions, 

J apibp2cp3dp4,ekifk2 \i=1..4 i=l,2 / 

X (/;(Pl)/;(P2)/;(P3)/^(P4)(l+/;(fel))(l+//(fe2)) - [/ O (1 + /)]) (2.15) 

^ Nevertheless the two previous references on bulk viscosity in a pion gas treated only elastic processes. 
Ref. [6] got a finite answer by using the methodology developed in [22], which assumes particle number 
is conserved and therefore allows a nonzero chemical potential in the equilibrium distribution function. 
Ref. [5] got a finite answer by doing a one-loop diagrammatic evaluation without "ladder" graphs, which 
amounts to neglecting the departure from equilibrium in all /'s other than /(p) in the Boltzmann equation. 
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The other contribution, C^"^^, arises when p = ki is one of the two pions. It is the same 
except 1^ is replaced with — ^jY?, so it cancels half of the above contribution. (These prefac- 
tors are symmetry factors to eliminate overcounting; for instance, if b, c, d are identical then 
only 1/3! of the phase space should be integrated over; and if 6, c, d are all distinct then the 
sum Ylibcd overcounts the possibihties by a factor of 3!. The sign difference arises from the 
relative sign between gain and loss terms.) 
Next we expand /^(p) to first order in ji, ST: 

W ^ Mp) + /o(p)(l+/o(p)) + ^) . (2.16) 

Inserting in Eq. (2.15) and expanding to first order in /x, 5T, the distribution functions 
become 

/o(Pi)/o(P2)/o(P3)/o(P4)(l+/o(fci))(l+/o(fc2)) (^(4 - 2)^ + (j^ - E^^) ^) ■ 

(2.17) 

The sum of energies cancels by energy conservation, leaving 

J apibp2cp3dp4,ekifk2 \l=l.A 1=1,2 / 

X (/o(Pl)/o(P2)/o(P3)/o(P4)(l+/o(fcl))(l+/o(fe2))) , (2.18) 

which determines /x. The value of n in turn determines the pressure correction, 

P-Pe^ = £ ^/o(p)(l+/o(p)) + ^) . (2.19) 

Recall that ST is set by the condition that the perturbation carry no net energy, which using 
Eq. (2.16) is 

ST /o(l+/o) = / (2.20) 

Jap Jap 

Together with Eq. (2.10) means 

P-Pe^^^iJ /o(l+/o)2 ^ 3^ . (2.21) 
Putting everything together with the definition Eq. (2.1), we find 



T{lapm+fo)2 



C = ^ , (2.22) 

4Cinel 



a 



inel 



4!2! 



\i=1..4 1=1,2 / 

X (/o(Pi)/o(P2)/o(P3)/o(P4)(l+/o(fci))(l+/o(fe2))) . (2.23) 



api bp2cp3dp4 ,efei/fe2 



The integration in the numerator is elementary, so evaluating the denominator will be our 
main challenge. 

Using the technique developed in [4, 10, 23], we would arrive at the same result by using 
the single parameter Ansatzioi the departure from equilibrium shown in Eq. (2.16). In the 
notation used there, each term in the numerator is S and the denominator is C. The factor 
of 4 is essentially {/i + ii + ii + ii — /i — / /i^ and can be understood as follows; each number 
changing collision changes particle number by 2 (one factor of 2), and a chemical potential 
makes the forward process faster than the backwards process by 2/1 /T (the other factor of 
2). 

III. CHIRAL PERTURBATION THEORY 

Quantum Chromodynamics is considered as the fundamental theory for describing strong 
interactions between quarks and gluons. However, at energies below the breaking scale of 
chiral symmetry, quarks and gluons are confined within the asymptotic hadron states, such 
as pions, kaons, and rj mesons. In this energy regime, the QCD coupling constant becomes 
so large that the theory is highly non-perturpative and we still lack an analytical method to 
solve it. However, the situation gets better if we write an effective field theory describing the 
meson states. It is an experimental fact that, at sufficiently low energies, the light mesons 
interact weakly with each other, with the strength of interactions controlled by a derivative 
expansion which is described by Chiral Perturbation Theory [18, 24], an effective theory for 
the interactions of light pscudoscalar mesons. 

In the chiral limit, the QCD Lagrangian possesses a SU {N)j^ x SU (iV)^ x U {l)y global 
symmetry. Here A'^ denotes the number of fiavors. The axial symmetry U (1)^ of the 
QCD Lagrangian, present at the classical level, is broken due to a quantum anomaly. Ex- 
perimental facts, such as the hadron spectrum and quark condensate, indicate that the 
SU (N)^ X SU (iV)^ X U {l)y spontaneously breaks down into SU {N)y x U {l)y. Accord- 
ing to Goldstone's Theorem, in this process, massless Goldstone bosons, which arc identified 
with the pseudoscalar mesons, are generated. Since we are dealing with a pure pion gas, we 
only focus on the case that N — 2, that is, only up and down quarks are of concern in our 
discussion. 

In this specific case, the three kinds of pions are considered as the Goldstone bosons, 
and they transform as a triplet under the subgroup SUv{2). Moreover, pion fields, the 
three-component vector $ = (0i, 02,03), are isomorphic to the quotient group SU (2)^ x 
SU{2)^/SUi2)y. 

In the chiral limit, one can, in terms of pion fields $ — (0i, 02, ^s), construct the general 
Langrangian invariant under SU (2)^ x SU (2)^ x U (l)y, with the ground state invariant 
only under subgroup SU {2)y x U (l)y. But in fact, instead of being massless, pions have 
small but finite masses around 135 MeV. This is because chiral symmetry is not an exact 
one. It is broken by a small amount due to the nonvanishing masses of up and down quarks. 
In order to give masses to pions, one also needs to add an explicit symmetry breaking term 
into the Lagrangian, which is treated as a small perturbation. 

The general effective Lagrangian can be organized by chiral order, 

J^eG — jO.2 + jO,4 + jCq + ■ ■ ■ 



7 




FIG. 1. Three classes of diagrams needed to evaluate the inelastic scattering rate to lowest order 
in chiral perturbation theory. Here the Roman subscripts are the Cartesian isospin indices and the 
number "2" in a circle denotes the chiral order of the vertex. 



where the subscripts indicate the chiral order. £2 with the smallest chiral order contains 
the minimum number of derivatives and quark mass terms. It reads [18] 

£2 = -^Tr {d^Ud'^U^) + lo^Tr {U + U^) . (3.1) 

Here 



T ■ $ 

U = exp ( i I = exp 



'.<Hx) 



(3.2) 



where Fq ^ 93MeV is the pion decay constant and r are the three Pauli matrices. 

The matrix element for elastic scattering is well known in chiral perturbation theory [18] 
and does not concern us, since Eq. (2.22) shows that the bulk viscosity is controlled by 
number-changing processes. We need the matrix element for Att — > 27r processes. Three 
classes of diagrams can arise, as depicted in Figure 1. For each class, we must sum over the 
distinct permutations of the external lines. 

Expanding £2, one can find the corresponding matrix elements. For the representative 
permutations shown in Figure 1, the matrix elements read (here for simplicity of writing 
down the matrix elements, all the four-momenta are viewed as incoming) 

Mi= V{a,b,e,g) ^ V{c,d,f,g) (3.4) 

M2= V{a,b,c,g) ~' ^ V{d,e,f,g) (3.5) 

■^3 = ^S^'S^'S'^ [4 {pa -Pb + Pc-Pd + Pe-Pf)- 3ml] 

+ all distinct pairings of the set {a, b, c, d, e, /} , (3.6) 

where Yl is a sum over the species type in the propagator, pg is the four-momentum of 

9=1,2,3 
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the propagator and 



V (a, /3, 7, g) = (^/3F2) [ 5'^^^^ (2p, • pp + 2p, • - % • + 

+ ^^''^^"'^ (2p„ • + 2p^ • p^ - 4p„ . p^ + m^) ] . (3.7) 
Therefore, the transition amphtude of the lowest order in question is 



\M\ 



(3.8) 



perm perm 

where ^ means a sum is taken over all distinct permutations of the external lines. 

perm 

This transition amplitude has a very complicated form, so we cannot finish the integral 
Cinei analytically. Therefore, we resort to numerical methods. For the numerical calculation, 
we work in the local plasma rest frame, that is, the rest-frame four velocity is — (1, 0, 0, 0). 
The distribution function in this frame is just /o (p) = [exp (Ep) — The main challenge 
is to perform the phase space integration over 6 external states. Wc consider the process 
as 47r — )■ 2tt, that is four incoming paticles and two outgoing particles. We perform uncon- 
strained integrations over the four incoming particle momenta in spherical coordinates with 
Pa as the z axis and pb lying in the x, z plane. 



d^Pad^Pbd^Pdd^Pe 1 f pldpaPldpbdeb'pldpbdVt^pldpbdVtd 



{2Try' 16EaEbE,Ed 8(27r)'V Ea E^ E^ E, 



(3.9) 



and then apply the energy-momentum conserving delta function to simplify the two-particle 
final phase space integration in the manner shown in [25, 26]. The final state phase space 
can be rewitten as 



/ 



^'P'^'P' + + + _ p. _ . v^l-y/' I (3.10) 



{27ry4EeEf 



where Q* is defined in the center of mass frame of the total incoming momentum k'^ — 
Pa+Pb ~^Pc ~^Pd^ and s — —k'^ is the Mandelstam variable. In the center of mass frame it is 
most convenient to work in spherical coordinates with the z axis chosen along the boost axis 
to the plasma rest frame. All dot products between incoming momenta are easily expressed 
in terms of the plasma frame variables, as is the Mandelstam variable s. For final state 
particle energies and dot products between an incoming and an outgoing momentum, we 
need to apply the boost between center of mass and plasma frame variables. An alternative 
approach is to consider the process 27r — > Air and apply the energy-momentum conserving 
delta function on the 4-particle final state phase space as shown in [25, 26]; but this approach 
is a little more involved. The resulting 11-dimensional integrations are performed by Monte- 
Carlo integration using CUBA [27]. 

We determine the pressure, speed of sound, and numerator of Eq. (2.22) by performing 
the integrals in Eq. (2.8), Eq. (2.9), and Eq. (2.22) numerically. It has become customary 



9 



T (MeV) 
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20 


30 


40 


50 
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70 




C (GeV^) 


3.6 X 10^^ 


2.1 X 10^ 


9.3 X 10^ 


3.9 X 10^ 


4.1 X 10° 


7.0 X 10" 


-1 


1.6 X 10" 


-1 


s (GeV^) 


2.4 X 10-^° 


3.9 X 10-^ 


5.8 X 10-^ 


2.7 X 10-^ 


7.4 X 10-^ 


1.6 X 10" 


-4 


2.9 X 10" 


-4 


T (MeV) 


80 


90 


100 


110 


120 


130 




140 




C (GeV^) 


4.7 X 10-2 


1.6 X 10^2 


5.9 X 10^3 


2.4 X 10-3 


1.1 X 10^3 


5.2 X 10~ 


-4 


2.6 X 10' 


'4 


s (GeV^) 


4.7 X 10""^ 


7.1 X 10-"^ 


1.0 X 10-3 


1.4 X 10-3 


1.9 X 10-3 


2.5 X 10" 


-3 


3.2 X 10" 


-3 



TABLE I. Values of C and s at certain temperatures 




20 40 60 80 100 120 140 j/^ 
T (MeV) 

FIG. 2. The numerical calculation of bulk viscosity ( and the bulk viscosity to entropy ratio (/s. 

to compare viscosities with the entropy density s — dP/dT, which has the same units as 
Differentiating Eq. (2.8), 

« = £^/o(l+/o) (3.11) 

which we also handle numerically. 



IV. RESULTS AND DISCUSSION 

The results of numerical calculation of the bulk viscosity are shown in the Table I and 
Figure 2. The most obvious feature of the bulk viscosity is that C. and C/s both rise as the 
temperature is lowered. This is the same behavior as the shear viscosity, in contrast to the 
high temperature regime, T ^ Aqcdj where r]/s rises but (/s falls with rising temperature. 

We can understand the rising behavior of (/s with lower temperature, for T ~ 111^^, as 
follows. First, the strength of conformal symmetry breaking depends on m-^/T, which gets 
larger as T gets smaller. Second, as the temperature gets lower, the typical momentum scale 
for pions gets lower. Since pions are pseudoGoldstone bosons, they interact mostly through 
high-derivative interactions, which get weaker as the energy scale is lowered. Therefore 
the system remains out of equilibrium longer, leading to higher viscosities. This last effect 
becomes very important when T <^ m.,r- In this case the density of pions falls exponentially, 
n ~ (m,r2^)^^^ exp(— mTr/T). The probability to have four pions in one place at one time. 
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to participate in a number-changing collision, is therefore exponentially small, ^ so the rate 
of number changing processes is exponentially suppressed and the bulk viscosity becomes 
exponentially large. This behavior was pointed out in the context of scalar field theory by 
Jeon [16]. 

In the low temperature limit T <^ m,r, the behavior of the bulk viscosity can be calculated 
analytically. In this regime the distribution function for incoming pions is well approximated 
by the nonrelativistic form fo{p) — e~"^/^e~^ /2m^T r^^ie typical value of the momentum p 
is p ~ \fm^ <^ m^r, which greatly simplifies both the initial particle phase space and the 
matrix element. For the purposes of evaluating the matrix element A^4_j.2, at leading order 
we can make the approximation that 

Po = Pb = Pc = = 0^ , Pe = (2m, v^m^ and p/ = (2m, -v^rn^ . (4.1) 

Under this approximation the summation of matrix element over species can be found in 
closed form: ^ 1-^1^ ~ 2025m^/2FQ^. Factoring it out of the integral, and approximat- 

ing s ~ 16m^, the remaining angular integrations can be performed easily. Then the phase 
space integral in Qnei reduces to 

X/o(Pa)/o(P6)/o(Pc)/o(Pd)(l+/o(Pe))(l+/o(p/)) (4.2) 

(4.3) 



40967r9m4 
163847r7 



We also need to carry out the integral in Eq. (2.22), which includes determining the speed 
of sound from Eq. (2.10). Here there is a subtlety; if we compute to lowest order and put 
it in Eq. (2.22), again computing in the nonrelativistic approximation, we get zero. Both 
equations must be expanded to second order in T/m^, yielding 



J^^/o(l+/o)2^1^^ = exp(-m./r)x (-3^|^ + 0(r'/X-''')) . (4.5) 



where the factor of 3 counts the number of pion species. Combining these results, the low 
temperature limit of the bulk viscosity is 

, 16384^377^ Fq^ 2m^ C,^ ^ 32768A/67r^Fo« 3m^ 

C(r«m)~ — |exp— ^, ^(r«m)~ ^^^exp— ^, (4.6) 

225 m^ i s ^ 



^ Or for the inverse process, the probability to have two pions with enough energy to generate four pions is 
exponentially small 
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— 1 3 

where we used the leading order behavior of Eq. (3.11), s ~ (3m^T2/(27r) 2) exp(— m^/T). 
These low temperature asymptotics are consistent with our numerical results. 

We should emphasize that at temperatures such that the bulk viscosity is very large, 
CV ■ u > P, the near-equilibrium expansion implicit in defining and using ( has broken 
down. When this occurs, the system in question has fallen out of chemical equilibrium; in 
fact CV • u > P can be taken as a criterion for the breakdown of chemical equilibrium and 
the freezing out of number changing processes. And when C becomes exponentially large, 
the approximation that we treat QCD without including electromagnetic interactions ceases 
to be valid. At low temperatures the dominant number changing process would actually be 
tt'^ — )■ 27 (and its crossings). We will not consider this extension here. 

In conclusion, we have computed the bulk viscosity of a pion gas, the natural low- 
temperature limit of QCD. We find that the bulk viscosity rises at low temperatures, growing 
exponentially as C ~ exp(2TO^/T) in the T <^ m^r limit. This growth implies that kinetic 
theory will generally break down at low temperatures, explaining chemical freezeout. 
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